Dynamics of precipitation pattern formation at geothermal hot springs 



O 

o 

<N 



<N 

in 

Oh 

d 



> 
oo 

O 

O 

o 



> 

• i-H 

x 



Nigel Goldenfeld, Pak Yuen Chan and John Veysey 
Department of Physics, University of Illinois at Urbana- Champaign, 
Loomis Laboratory of Physics, 1110 West Green Street, Urbana, Illinois, 61801-3080. 

We formulate and model the dynamics of spatial patterns arising during the precipitation of 
calcium carbonate from a supersaturated shallow water flow. The model describes the formation of 
travertine deposits at geothermal hot springs and rimstone dams of calcite in caves. We find explicit 
solutions for travertine domes at low flow rates, identify the linear instabilities which generate dam 
and pond formation on sloped substrates, and present simulations of statistical landscape evolution. 
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The terraced architecture of carbonate mineral de- 
posits at geothermal hot springs is one of the most strik- 
ing and beautiful terrestrial landscapes. Spring water at 
above 70° C emerges from the ground at a vent, releases 
CO2 and precipitates CaC03 in the form of travertine, 
as it flows downhill over the pre-existing terrain [J, 0, 0. 
The terrain itself is thus constantly changing in response 
to the influx of CaC03, with measurements indicatin 
precipitation rates as high as 1-5 mm per day0, 0, 
The ever-changing substrate modifies the flow path of the 
spring water, resulting in a constant dynamic interplay 
between the landscape and the fluid flow as the travertine 
outcrop develops and grows. The resulting morphology 
is a cascade of nested ponds and terraces at a wide va- 
riety of scales ranging from hundreds of meters down to 
millimeters, as shown in Fig. Q}. 

Most studies of this phenomenon have tended to fo- 
cus on the microscopic origins of the crystal growth pro- 
cess: what is the extent of biomineralization due to ther- 
mophilic microbes 1^, what are the primary controls 
on degassing || 0, mineral composition [1(1 Hl| . crystal 
fabric and habit |l2^. and structure [13. Il4j? These issues 
are complex and system dependent, but a main theme is 
the competition between biotic and abiotic mechanisms 
of precipitation. Microbial metabolic activity can in prin- 
ciple locally influence the CO2 composition of the water 
column, thereby affecting the rate of CaCOs precipita- 
tion. However, turbulent flow over pond lips or evap- 




FIG. 1: (Color online) Travertine formation at Angel Terrace, 
Mammoth Hot Springs, WY. (a) a large pond, of order 1 
meter in diameter, and smaller features, (b) a portion of the 
flow system about 25 meters from the vent, on the scale of 
centimeters. 



oration can be more effective degassing mechanisms 9j. 
In fact, the issue of biotic versus abiotic mechanism is 
not straightforward, because even if the main mechanism 
for degassing is abiotic, the kinetics of precipitation and 
crystal growth may require the presence of exogenous 
particles (such as dead or living bacteria) to initiate nu- 
cleation. 

The purpose of this Letter is to address instead the 
origin of the large scale structure of the terraced archi- 
tecture. The ubiquity of this morphology at carbonate 
geothermal hot springs world-wide, as well as at low tem- 
perature speleothem rimstone dam formations 15], sug- 
gests that it is appropriate to seek a generic explanation 
based on principles of fluid dynamics, precipitation kinet- 
ics, and crystal growth dynamics, rather than one that 
hinges on specific material parameters or system com- 
ponents such as microbes. In fact, we will see that it is 
possible to explain in this way the formation of ponds, the 
variations in terrace morphology, the absence of an obvi- 
ous characteristic lengthscale of the landscape, and even 
the quantitative properties of simple landscape motifs, 
such as the circularly-symmetric travertine dome shown 
in Fig. (Hi). 

The formation of a terraced architecture is a result of 
a depositional instability arising from turbulent flow of 
a supersaturated solution over a surface, and represents 
an example of free boundary dynamics in precipitative 
pattern formation, related to the phenomena that give 
rise to stalactites |ld Il7j. 

Travertine domes:- We begin by analyzing the dynamics 
of pattern motifs, ignoring interactions between the m, by 
analogy with earlier work on solidification patterns 1 181 . 
In the present problem, these features are pond lips|l£j 
and travertine domes. 

We now formulate a boundary-layer model of the 
growth of travertine domes, coupling the evolution of the 
travertine substrate to the fluid dynamics in a thin film 
of depth h around it. The kinematic equation govern- 
ing time evolution of the curvature, ft, of a curve in two 
dimensions is given bv[I3. Eo|: 
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FIG. 2: (Color online) Travertine dome at Mammoth Hot 
Springs, WY. (a) Dome whose central pond is 50.3cm in di- 
ameter, (b) Dome height as a function of 0. Horizontal axis: 
analytical prediction from Eq. 0). Vertical axis: simulated 
dome (CA) with no surface tension (squares), with surface 
tension (circles), and the profile of the dome in (a) (trian- 
gles, data digitized from photograph). The points all fall on a 
straight line with unit slope above the fluting point, indicat- 
ing agreement with Eq. fitting only ro. Below that one 
extra parameter is required, as explained in the text. Inset: 
comparison with dome from (a). 



where 6 is the angle between the surface and the ver- 
tical axis and v n is the normal growth velocity of the 
surface. The time derivative in the equation is defined 
with respect to fixed 0. This equation is purely geo- 
metrical; for any given function v n of water chemistry, 
surface kinetics, fluid flow state, the evolution of k is de- 
termined. Here, we follow Wooding ^| and make the 
simple assumption that v n is directly proportional to the 
depth- averaged tangential fluid velocity U: v n = GU 
where G is a mass transfer coefficient, depending on the 
water chemistry and the turbulent flow near the growing 
surface jjl). 

In general, we need to couple the Navier-Stokes equa- 
tion to Eq. JU to obtain a complete description of 
the coupled fluid dynamics and surface kinematics, but 
there are a number of simplifications. First, because the 
growth rate is of order 1 — 5mm/day and the fluid flow 
rate is of order lmm/sec, there is a separation of time 
scales. So if we are interested in the morphological evo- 
lution, we can neglect the change in flow rate, represented 
by the time derivative in the Navier-Stokes equation. 
Second, our field observations indicate that the thickness 
of the fluid film flowing over the domes is very small com- 
pared to the curvature of the surface; thus, we make the 
approximation that the fluid is flowing down a (locally) 
constant slope. In addition, the flow is apparently lam- 
inar, so that we can use the Poiseuille-Hagen profile for 
the velocity in thin film to give a depth- averaged mean 
flow velocity U : 
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where a = gQ 2 /12ttu : g is the gravitational acceleration, 
Q is the total mass flux coming out of the vent, v is the 
viscosity of the fluid and r is the axial distance from the 



vent. Circular symmetry is imposed to arrive at J2J). We 
will later see that the assumption of laminar flow is self- 
consistently verified. For a dome, steadily translating 
upwards without change of shape with velocity v t , Eq. 
(PJ gives 
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Rearranging terms gives the shape of the dome as a one- 
parameter family of curves 



r(0)/r o = 



sin# 
cos 3 
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where the scale factor ro = \/G 3 a/vf. Eq. (@J) is 
plotted in Fig. (JSb). Good agreement is obtained be- 
tween our theory and the observations below a critical 
angle C . From the fit, and the typical parameter values 
G ~ 10 -8 , vt ~ lmm/day and Q ~ lcm 3 /sec, we obtain 
U ~ 25mm/sec and h ~ 1 — 10mm, and a Reynold's num- 
ber, Re = Uh/v ~ 10 — 100. The assumption of laminar 
flow is self-consistently verified. 

For angles > C , the analytical profile does not fit the 
field observations; the point of departure closely follows 
the point where we also observe a fluting pattern around 
the dome. We will show below that this is due to the 
effects of surface tension at the air-water-travertine in- 
terface. As the water flows out of the dome, it is spread 
over an increasingly larger area and thus the fluid thick- 
ness decreases, ultimately reaching a point where contact 
lines form and surface tension can not be ignored. The 
analytical solution above neglects surface tension, but 
is able to predict the limiting point of its validity when 
h ~ d Cl leading to a prediction for the scaling dependence 
of the critical angle on the model parameters. 

The inclusion of surface tension introduces an addi- 
tional length scale, namely, the capillary length, d Cl into 
the problem. Now, the only other length scale in the 
problem is ro = y gG 3 Q 2 jvv\. Since C is dimensionless, 
it can only depend on the ratio ro/d c and G. For a given 
chemical environment, G is fixed and we are left with the 
prediction, derived from our analytical solution, that 
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verified below in Fig. (03). 

The damming instability:- To understand the lack of any 
characteristic scale in the landscape, we consider the sta- 
bility of the moving boundary problem for turbulent fluid 
flowing down a constant slope, on which deposition may 
occur. In order to capture the turbulent flow, we make 
two approximations. First, we use the thin film approxi- 
mation attributed to de St. Venant, which is valid when 
the fluid film thickness is much less than the character- 
istic scale of variation of the flow in the streamwise di- 
rection, but include the lowest order corrections for the 
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curvature of the underlying surface [U EH- 

d t u + d s E = -C f gu 2 /gh(l - kK) , , 

(1 - Kh)d t h - d s q = w 

where 

u(s,n,t) = uo(s,t)/(l — kJi), 

E(s,t) = C J rhco^0^ul/(2g(l-nh) 2 ), (7) 
q(s, t) = -u / (k, log(l - Kh)), 

where u is the fluid velocity, h is the fluid thickness, £ 
is the height of the underlying surface measured from a 
fixed horizontal axis and s is the arc length measured 
from the top of the system. Secondly, we have modeled 
turbulent flow phenomenology by including the term pro- 
portional to Cf, the Chezy coefficient |24j . which empir- 
ically describes the energy lost due to turbulence, in a 
manner consistent with Kolmogorov's 1941 scaling the- 
ory of turbulence (K41)|25L l26j |. These two equations 
have to be solved together with the growth equation (QJ. 
The trivial solution to this set of equations can be easily 
found, and is simply uniform viscous flow down a slope. 

To study the linear stability of this solution, we add a 
perturbation proportional to exjp(ikx + cu(k)t), and cal- 
culate the spectrum of the growth rate a; as a function 
of the wavenumber k for the linearised set of equations. 
It is found that Re(cj) is positive for all values of fc, indi- 
cating that the system is unstable towards perturbations 
of all length scales, with no special scale singled out at 
linear order. 

We next present an approach to simulate the statis- 
tical correlations of the terraced landscape in the fully- 
nonlinear regime. 

Cellular modeling:- We represent the landscape as com- 
posed of stacked "bricks" or cells by its height H(i,j), 
above a horizontal reference plane, where i and j are 
x — y coordinates in the reference plane. The water col- 
umn is situated above the height field and represented 
by the variable W(i,j) describing the volume of water 
above each coordinate element. Each packet of water 
can also contain calcium ions, C(i, j), and dissolved car- 
bon dioxide vapor, V(i,j), which may potentially cause 
precipitation through a caricature of the complex reac- 
tion pathway given by Ca 2+ + 2HCO3 ^ CaC0 3 (s) + 
H2O + C02(g). The evolution of the landscape is gov- 
erned by update rules on these fields, which mimic and 
go beyond the continuum description given above. In 
addition, at each point, it is necessary to keep track of 
water that is ponded, W p (i,j), and the temperature of 
the water T(i,j). 

A complete lattice update consists of the following 
steps, which we will describe in more detail below: (1) 
Add water to the system, at the spring source taken to 
be the origin; (2) Propagate all the water in the system, 
by moving packets to nearest and next nearest neighbour 
grid points, and ensuring that the water in all ponds is 



level; (3) Update the water chemistry, e.g. C(i,j) and 
V(i,j) to take into account outgassing due to fluid mo- 
tion and depletion of Ca ions due to precipitation; (4) 
Evolve the height field in response to the precipitation of 
CaC0 3 . 

In step (1), a quantity of water SW is added at the 
source, so that W f (0) = Wo, a constant value appropriate 
for a constant pressure head (neglecting the change in 
pressure due to the vertical growth of the landscape). 
Here and below, primed quantities denote the updated 
variables. The new water added to the system contains 
initial concentrations of calcium, Co, carbon dioxide, Vb, 
and is at a temperature To; the fact that the source water 
is undersaturated is represented by Co < Vq. The values 
of these fields at the source point are updated based on 
the volumetric ratio of the amount of water added to 
the existing water: e.g. C = (C(O)W(O) + 5WC )/W , 
and similarly for the CO2 concentration and temperature 
fields. 

The transport of water in step (2) is carried out by a 
variation of an algorithm used for braided river flow 27J, 
in which the flux along bonds connecting a given lat- 
tice point to one of its eight closest neighbours is de- 
termined by the landscape gradient along that direction, 
while conserving the total volume of water. If the slope 
S exceeds a threshold 5 C , the flow is taken to be turbu- 
lent and the flux is proportional to y/S in accord with 
Chezy's law. The appropriate height variable for deter- 
mining the chemical potential of the water, and hence 
the equilibrium filling of ponds in the landscape is the 
"energy surface" Ht = H + W. Water moving on this 
surface is also subject to surface tension and contact line 
effects, that can be important near the rim of a pond, 
for example [23. This is modeled by requiring that wa- 
ter is propagated only if W exceeds a small threshold 
at that point. Packets of water carry with them ad- 
vected variables T, C and V, which are updated by the 
volumetrically- weighted average of all the neighbourhood 
points from which the water packet originated. 

Water chemistry is updated in step (3), by allowing 
a constant fraction of CO2 to outgas at each time step, 
with an additional outgassing component proportional to 
y/~S / (l-\-y/S) to reflect the influence of slope-initiated tur- 
bulent flow. As CaC03 is deposited, the concentrations 
C and V change to reflect mass balance. 

The evolution of the height field H is the final step 
in the lattice update, with a change given by SH = 
(C - V) x (R ± + R 2 F • n + R 3 (S)). Here R x and R 2 
are positive constants, F is the flux between cells, with a 
direction given by the component of the gradient of the 
energy surface Ht between cells and magnitude given by 
the volume of water propagation per unit time step, ft is 
the unit vector normal to the underlying surface H, and 
Rs(S) is proportional to y/S, representing the increased 
precipitation due to Bernoulli effects and local turbulent 
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FIG. 3: (Color online) Results from simulation of the CDS 
model, (a) Portion of a typical simulated landscape. The fig- 
ure clearly shows distinct geomorphological regimes or facies, 
as observed in the field 6] . (b) Critical angle for the contact 
line formation on a travertine dome, plotted according to Eq. 
®, showing data collapse, as predicted by theory. Inset: raw 
data. 

degassing. 

In Fig. (J3t) is shown a snapshot from a typical time- 
dependent simulation, initiated on a sloping plane with 
small initial roughness. For generic values of the model 
parameters, we observe the depositional instability pre- 
dicted above, and the formation of ponds and terraces 
in broad qualitative agreement with field observations. 
The same calculation also yields travertine domes when 
initiated on an initially horizontal surface. Fig. (03) 
shows that the CA model agrees with the analytical the- 
ory when the surface tension is switched off. Moreover, 
the same CA model predicts the exact shape of the ob- 
served travertine domes when surface tension is switched 
on, with a contact line at C where fluting emerges. The 
scaling prediction Eq. J5J) is verified in Fig. J3t>) over 
5 decades of Q 2 v^/dl; the inset shows the raw (un- 
sealed) data for 6 C as a function of d c for v t ~ 0.5 and 
0.1 < Q < 5000. Our results show that travertine pre- 
cipitation pattern formation results from an interplay be- 
tween fluid flow, capillarity and chemistry, that can be 
captured quantitatively from both analytical and cell dy- 
namical system approaches. 

Future field work and extended simulation studies will 
explore the statistical properties of terraced landscapes, 
and whether these have a universal character computable 
with no adjustable parameters by minimal models. For 
example, layer fluctuations in ancient stromatolite rocks 
exhibit power-law correlations accurately captured by a 
simple generic model ji^ in the same spirit as the one 
given here. Our preliminary simulations indicate that at 
least some statistical measures of the landscapes, such as 
the pond size distribution function, are fractal. 
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